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ABSTRACT 

We determine approximate numerical integrals of motion of 2D symmetric Hamiltonian systems. We detail for a few gravitational 
potentials the conditions under which quasi-integrals are obtained as polynomial series. We show that each of these potentials has a 
wide range of regular orbits that are accurately modelled with a unique approximate integral of motion. 
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1. Introduction 

Gravitational potentials with explicitly known second integrals 
of motion are rare, the most freq uently used in astronomy are 
axisymmetric potentials and also Stackel ( 1893 i potentialsMthat 
cover a large class of potentials, but are not always sufficiently 
realistic. Such explicit forms are useful, for instance, to under- 
stand more clearly the underlying physics of dynamical systems, 
to build stationary distribution functions, etc. |Hietarinta| (1987 ) 
gives such a list of other potentials with known second integrals, 
but potentials with an analytic second integral of motion are the 
exception as shown by Poincare, and ergodic motions should be 
the rule. However, according to the KAM theorem, at low ener- 
gies many orbits are regular and remain confined on tori, while 
for irregular orbits the time diffusion may be exponentially long 
( |Nekhoroshev| 1 977| |Mo"rbidelli & Giorgillil 1 995a"|b"l >. For these 
orbits, as shown numerically by Ollongren ( 1962) for a realistic 



galactic potential, we expect the existence of approximate first 
integrals. 

The need for generality leads to numerous efforts to ob- 
tain tractable integrals of motion for more realistic potentials. 
Among the developed methods, we can mention the dynamical 



spectra and more specifically the frequency analysis (Binney & 
Spergel 1982; Lask ar|1993 i with a recent 3D application to nu 



merical simulations of galaxies ( |Valluri et al.|2012| . Frequencies 
associated to each orbit are constants of motion and allow accu- 
rate classifications of orbits. Another example, closer to the work 
developed here, is the tori reconstruction of M cGill & Binney| 
(1990), who modelled the tori on which orbits circulate in the 
phase space. Their analysis is based on the angle-action vari- 
ables, and the time dependency remains explicit. This method is 
very attractive since the dynamics of the system becomes very 
simple. They distort analytic tori of a toy Hamiltonian into the 
tori of the Hamiltonian of interest. This allows the analysis of 
any non-integrable potential as a perturbation of a nearly inte- 
grable one (Binney & Kumar 19931. Extensions for the sepa- 



1 Stackel potentials, as they are called in the astronomical literature, 
are usualy named Darboux potentials in the mathematics literature. A 



first systematic attempt to study quadratic invariants is given in D arboux| 
(1901). Contributions on separability, orthogonal coordinate systems, 



rate modelling of orbit families are described by |Kaasalainen| 
|& Binney| ( |1994| ), |Kaasalainen| ( |1994||1995[). More direc t mod- 
ellings of orbits are also proposed by [Prendergast ( 1982) using 
the ratio of trigonometric functions, or by |Robnik| ( |1993] ) us- 
ing Pade approximations (rational functions). Close to our work, 
Warnock (1991) determined approximate integrals by fitting or- 



bits. He developed a fit using an approximate discrete Fourier 
transform of positions along orbits and recovered the coefficients 
of an integral of motion. This method works for non-resonant or 
high-order resonant orbits but fails for low-order resonant or- 
bits. |Bazzani et al.| ( |1991| ) determined an approximate integral 
for (2D) sympletic maps, minimising the variation of the approx- 
imate integral along the orbit under study. 

Formal integrals of motion can be obtained by directly solv- 
ing the Boltzmann equation assuming that the integral may be 
written as a polynomial (see Contopoulos 1960| ). The conver- 
gence of these series is not guaranted, but they may be asymp- 
totic (i.e. semi convergent series), allowing one to approximate 
integrals. Birk hoff| ( |1927[ ) and Gustavson (1966 ), using normal 
forms, proposed a process to build such approximate formal in- 
tegrals of motions. A more complete description of all these sub- 
jects can be found in Contopoulos (2002). Introductory lectures 
on these questions of galactic dynamics and celestial mechanics 



may be found in Henon ( 1983) or Efthymiopoulos et al. ( |2007 



and quadratic invariants may be found in Ankiewicz & Pask 1 1983b. 



The motivation of the work presented in this paper is to ob- 
tain tractable and approximate integrals of motion for a few sim- 
ple potentials that are representative of galactic potentials. We 
restricted our study to 2D motions. Future applications might 
be the building of dynamically self-consistent galactic mod- 
els for stationary potentials, and modelling the kinematics of 
its stellar populations with distribution functions. These goals 
are usually achieved with N-body numerical simulations, with 
the |Schwarzschild ( 1979 ), or with the Syer & Tremaine 
methods. 

Our method consists of solving the stationary Boltzmann 
equation for a second integral of motion. Instead of using a for- 
mal integration, we numerically determine the coefficients of a 
series by a least-squares minimisation of the integral variance 
along orbits. This allows us to obtain a single expression for 
different orbit families. The method is described in Section 2 
and detailed explanations and discussion for a simple 2D po- 
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tential are given in Section 3. Applications to other potentials are 
briefly described in Sections 4-5, and the conclusion is presented 
in Section 6. 

2. Method and algorithm 

According to a definition proposed by Hietarinta ( 1987j), the 
concept of integrability is to be able to make some quantitative 
general statements of the dynamical system under study using a 
quantity whose value is conserved during the time evolution of 
a system. 

In this paper we approach this problem numerically and con- 
sider two degrees of freedom time independent Hamiltonian sys- 
tems of the form 

H= l -(u 2 +v 2 ) + (5(jc,y), 

with x, y the space coordinates and u, v their conjugate mo- 
menta. The energy is naturally a first invariant, and we search 
numerically for a second invariant I(x, v) = I(x, y, u, v) indepen- 
dent of the energy. 

We recall that this problem is closely related to the colli- 
sionless Bolztmann equation and that any distribution function 
independent of time of the form fix, v) = g(I(x, v), E(x, v)) will 
be the solution of the stationary 2D Boltzmann equation, 

df df d<&df d®df 

U h V = 0, 

dx dy dx du dy dv 

since this equation is also the Poisson bracket {/, H}=0 that is 
a property satisfied by any invariant of the Hamiltonian systems. 

From Hamilton's equations of motion, we can determine or- 
bits and note that if a second invariant I(x, v) independent of the 
energy exists, it remains constant, by definition, along any regu- 
lar orbit. Thus we assume the existence of such a second invari- 
ant, and we also assume that it may be written as a polynomial 
finite series of the coordinates and momenta (x,y, u, v): 

7 = /(*,v) = J] c Uv , n x k y l u m v\ (1) 

k,l,m,n 

To determine the coefficients of the series, we select a set of reg- 
ular orbits for a given potential, for each orbit a set of positions 
(x,y, u,v) along that orbit, and for each position I is evaluated. 
Coefficients are computed by a least-squares minimisation that 
minimises the variation of I around its mean value, (i.e. min- 
imising the standard error of /) along each orbit. This standard 
error would be zero if an exact integral I existed for the consid- 
ered potential and, according to our construction here, if it were 
also a polynomial. Before realising the minimisation, a number 
of considerations allows us to cancel redundant and useless coef- 
ficients within Eq. 1 . First, to avoid obtaining the trivial solution 
where all coefficients are null, we fix a coefficient to the value 
one: here, we will fix the coefficient of one these terms: either 
x 2 , v 2 , or x 2 v 2 . Many degrees of freedom remain to determine the 
coefficients within the series of I, and we may want that the pro- 
cess to build / leads to a unique solution. The unicity of adjusted 
solutions is necessary and useful to allow the comparison of so- 
lutions obtained using different orbits, different integration time 
for orbits, different number of coefficients, etc. We also want 
that the integral of motion I is independent of the energy. For 
that purpose, we modify I by removing the term u 2 '" from the 
series by subtracting from / a quantity proportional to E m (thus 
the modified / remains an integral). Starting with m — 1 and in- 
creasing m, we iteratively modify /, removing all terms u m , and 



now I(x, v) and E(x, v) are independent. All this is to explain that 
we do not remove useful terms, because the independence is sat- 
isfied if the rank of the Jacobian d(E, I)/d(xi, pj) is 2 (instead of 
removing the u m terms, we could have cancelled the coefficients 
of the v 2 " terms, or as efficiently those of the x k oiy 1 terms.) 

Some other coefficients are redundant. If / is an invariant, 
so is g(I), where g is any analytic function. To ensure that our 
minimisation process produces a unique solution, we reiterate 
and now subtract from / a quantity proportional to 7* (k > 1) 
to remove the term x k . Iterating this operatiorj^J all terms of the 
form x k , except the first one, are removed from the series. If the 
potential is even for the x coordinate, the first remaining term is 
in general the x 2 term. We find that selections of the coefficients 
and energy independence are critical for obtaining a numerically 
accurate second integral of motion. 

Finally, a significant number of other coefficients are re- 
moved for symmetry reasons. Each of the x and y axial symme- 
try of the potentials considered in this paper allows us to cancel 
half of the coefficients. Moreover, it is also known (Hietarinta 
1987) that if 7 is an integral of motion, it can be split according 
to the momenta parity, and the odd and even part 7+ and /_ are 
also integrals. If the potential is not superintegrable (and thus 
has at most two integrals), 1+ or /_ must be equal to zero (oth- 
erwise they are dependent), and then the momenta parity of / is 
clearly defined (in practice if the potential contains box orbits, 
the momenta parity of 7 is even.) 

Practically, we rewrite the polynomial 7, where the a, are N 
distinct monomials and the c,- are N adjusted coefficients, as 

I — flo + ^i=l c i a i ■ 

We set, for instance, ao = x 2 and its coefficient co is 1. 
From M positions along a given orbit, we set a im = 
di(x m ,y m , u,„, v m ), with m — 1 to M. 

We define the mean value of 7 over the M positions as 

and its standard deviation as 

Minimising the standard deviation reduces to N linear equa- 
tions equivalent to the matrix equation 

D.C + B = 

with 

Ci = ci , 

and 

— r J_y M l 

In the case of simultaneous fitting of P different orbits with 
M positions along each orbit, the linear equations remain similar. 
If the index m and p refer to the position m on orbit p, we have 

Bi — ^p=l Bi,p — ^ p=\ a 0,m,p - a i,m,p > 

2 Here, we do not answer the question of the convergence of these 
iterative operations. Our success in obtaining approximate integrals in 
the next paragraphs indicates that the iterations can build at least asymp- 
totic series. 
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and 



at 



M 



®i,m" ,p] • 



When M x P is larger than N, the problem is overdeter- 
mined, and we use the dgesv routine (LU decomposition) from 
the LAPACK software (www.netlib.org/lapack/) to solve the sys- 
tem of linear equations. Other LAPACK routines have been tried 
and give similar accuracies. 

The various fixed or removed coefficients depend on the 
shape of the potential. In Sections 3, 4.1, and 4.2 we set the fol- 
lowing coefficients: 

C2,0,0,0 = 1 

C2*,o,o,o = 1 with k > 1 
co,2/,o,o = with I > 1 



CkU 



if I + m or k + n or m + n is odd. 



In Section 4.3 the following coefficients are fixed: 

C0,0,0,2 = 1 

co,o,o,2« = 1 with n > 1 
co,o,2m,o = with m > 1 



CkU 



if / + m or k + n or m + n is odd. 



In Section 5 the following coefficients are fixed: 

c 2,0,0,0 = 1 

co,o,o,«2 = with n > 1 
C2*,o,o,0 = 1 with k > 1 
Ck,i,m,n = if k + n is odd orm + n is even. 

In summary, the method presented here to solve the 
Boltzmann equation and to obtain a first integral I(x,v) con- 
sists of computing some regular orbits, and of adjusting I(x,v) 
to these peculiar solutions. The constraint on the coefficients of 
the polynomial series modelling I(x, v) is that / remains constant 
along each selected orbit. The most important details are given 
within the next section using an 'exponential' potential. Other 
results are summarised for a few potentials in the following sec- 
tions. 



3. Exponential potential 

We apply the method described in Section 2 to obtain an approx- 
imate integral of motion to the case of the following barred po- 
tential: 

<&(x,y) = - exp(- q 2 x 1 - y 2 ). 

This potential is chosen for its simplicity: bounded orbits 
have energies in a limited range within [-1,0] and most orbits 
are regular. The density associated to the potential is not positive 
everywhere but the Taylor series of this potential has an infinite 
radius of convergence that should allow us, in a future work, an 
easier comparison of our results to formal integrals obtained by 
other means (Giorgilli & Galgani 1978 ). 

Most of the bounded orbits (E < 0) are regular. They are box 
and tube orbits (libration and circulation), and most of the tube 
orbits are of type (1:1). This can be seen in Poincare surfaces 
of section (Figs. [T| or with the 2D isolevel plot of x-apocentre 
(at u = y = 0) versus the energy and the initial position Xi n ; t (at 
u = y — 0) (Fig. [2l. In this last figure, high-order resonances 
appear at the highest energies and draw inclined 'plumes'. 

To build the polynomial quasi integral of motion, we select 
orbits that cover a wide range of energies and a wide range of 
second constants of motion, taken as the orbit initial position 




Fig. 1. Poincare surfaces of section for the "exponential" poten- 
tial with q = 0.8 and energies E = -0.5 (top) and -0.2 (down). 



*init (u = y = 0) (in this example, we consider only orbits that 
cross the x-axis perpendicularly, which excludes very few orbits 
that appear mainly at high energies.) 

We also select the polynomial terms used to determine /. Due 
to the x- and y-axis symmetries of the potential, the coefficients 
of terms in Eq. 1 with I + m or k + n odd are set to zero. The 
momenta parity of / is even (m+n even). We set the coefficient of 
x 2 to one, remove all x 2k (k > 1) and all v 2 ' terms. The final series 
thus obtained will be independent of E and a unique solution for 
the coefficients is expected. 

The remaining coefficients of the series in Eq. 1 are obtained 
by a least-squares minimisation, and the quality of the fit is quan- 
titatively defined by the constancy of / along each orbit. The re- 
sult of the fit critically depends on the existence or non-existence 
for the examined potential of an approximate integral that re- 
mains nearly constant over long periods of time along orbits. 
But the quality, also critically, depends as any polynomial fit, 
on a convenient coverage of the fitted space to adjust structures, 
and on a sufficiently large number of coefficients to model the 
numerous and small structures. Finally, it also depends on a suf- 
ficiently large number of fitted data to avoid or at least to min- 
imise erratic oscillations between fitted data. For this, we impose 
that the number of fitted positions is sufficiently large by apply- 
ing the recommended conditions that 2 » A^oeff ( |Dalquist| 
& Bjork 1974). We must also ensure a sufficient coverage: 1000 
positions on each orbit corresponds to a mean distance of 1 1 deg 
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Fig. 2. Isolevels of x-axis orbit apocentres versus Xmt an d E. At 
low energies box orbits dominate the diagram. Tube orbits (type 
1:1) become significantly present at E > -0.4. At E — -0.05, 
they cover the range Xi n j t from to 1 .4. Higher order resonances 
draw visible "plumes". 



iiiigiii 



Fig. 3. Second integral of motion / (polynomial of order 12) 
versus x; n i t for 2085 orbits covering 19 energies. Orbits with the 
same energy draw 19 dotted lines. Box orbits have / > 0, tube 
orbits have / < 0. The red continuous line is the dispersion of /, 
it increases at resonances. 



for the angle positions, assuming the couple of angles and ac- 
tions were known. 

Orbits are computed with a fourth-order Runge-Kutta 
integrator, the time step is fixed and the energy is conserved 
at 10~ 8 . The minimisation is written as a linear least-squares 
fit, inversion is performed using the LAPACK softwares and 
presents no peculiar difficulties, because the matrix size remains 
small. 



3.1. Results 

The first example is obtained with the following conditions: we 
select 2085 orbits covering energies E - -0.95, -0.9, -0.05 
and at each energy with initial conditions (jCinit, yinit=0, «init=0, 
Vi n i t (i?, Xinit)) Xi n it covering the interval to x max (E) by step 




Fig. 4. For 19 orbits with E = -0.05 close to the resonance 2:3, 
I{t) is plotted every At = 4 over the time interval T = 4000. The 
dispersion of I{t) increases for orbits close to the resonance. The 
time sequence is shifted for each orbit for visibility. 



Ax^ = 0.01. For each orbit, 1000 positions are taken (time step 
AT = 4) from T = to 4000. After removing cancelled terms, 
we use the first 200 coefficients of Eq. 1, a 12' /! order polyno- 
mial, and after fitting, following the procedure described in the 
previous section, we obtain /12 a quasi invariant integral of mo- 
tion. 

We find that the mean value, varies from orbit to orbit 

between 1 and — 1-4.6 (Pig. [3]» and that the mode (histogram 

maximum) of the dispersions cr /|2 along each orbit is 0.002, and 
<Tj remains within [5.10~ 4 - 3.10" 3 ] for the low-energy orbits 
(E < -0.5). At higher energies, 30 orbits close to resonances 
have 07 within [0.03 - 0.15], the remaining 1360 orbits have 07 
smaller than 0.03 and a mode equal to 0.002. 

We find that I\2 remains nearly invariant over long periods by 
computing for each of the 2085 orbits I\2 over the time interval 
AT = 4.10 6 , corresponding approximately to 50000 rotations, a 
much longer time than the interval AT = 4.10 3 used for the fit 
(50 rotations). The variations of /12 (i.e. its residuals) are identi- 
cal to a few percent to residuals obtained within the fitted time 
interval and, thus, we conclude that I\2 can be extrapolated in 
time far outside the fitted domain (but of course not outside the 
fitted domain in coordinates and momenta). 

We extend this check to examine the variation of In for 
orbits not used for the fit. Therefore, we explore in detail the 
energy-x,„„ fitted domain with a finer grid in energy (AE = 0.01 
and Axi n i t = 0.0025) using ~32000 orbits. For these orbits the 
dispersion remains very close to the dispersion of the nearby or- 
bits used for the fit. This confirms that we use a sufficiently large 
number of orbits and positions and that the interpolation between 
fitted orbits is correctly achieved. 

Figure[3]presents, for each of the 2085 orbits, the mean value 
of I\2 along the orbit versus its initial positions Xmt- For each 
orbit, Xj n i t is either its x-pericentre or its x-apocentre along the x- 
axis (thus with y-u-Q). Each dotted line represents a sequence 
of orbits with an identical value of the energy. 

The two main families of orbits are box orbits and tube orbits 
(with type 1:1), corresponding to positive and negative values of 
I\2- The family of box orbits is dominant at low energies, E S> 
-0.4. At all energies, the 1:1 periodic tube orbit corresponds to 
the minimum of /12. For most of the orbits, there is a one-to- 
one relation between each orbit and a couple of values (E, I). 
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Fig. 5. - Exponential potential, (x, u) Poincare surfaces of section at E — -0.05 from the numerical and analytical /: (top-left) 
Numerical section, (top-right) analytical section obtained from a fit at order 12 and orbits at energies from -1 to -0.05, (bottom-left) 
from a fit at order 18 and orbits with energy E = -0.05, (bottom-right) from three local fits (order 18) of orbits in the neighbourhood 
of three resonant orbits, the empty space between the islands of stability is not filled correctly by the analytical integrals (red dots 
are the sections, y — and v > 0, from computed orbits). 



This useful property is lost, however, when the fit is improved to 
adjust higher order resonances, as we show below. 

As a consequence of our specific choice of the cancelled co- 
efficients within Eq. 1, we have / = for the radial orbits aligned 
along the y-axis. On the other side, the x-axis radial orbits are, 
at a given energy, the orbit with the largest x-apocentre, X sup ^, 
and we have I(X sup , E ,0,0,0)=X 2 mp E . 

In Figure [3] the dispersion <x /p is also plotted. This disper- 
sion increases abruptly at the proximity of extended resonances. 
In Figure|4j we plot /12(f) for 19 orbits close to a resonance 2:3 
at E — -0.05. In this figure, for each orbit /12(f) is plotted every 
AT - 4 from T = to 4000 (and each plotted orbit is shifted 
for better visibility.) Crossing the resonance, the dispersion in- 
creases and is of the order of the variation of /, and within the 
resonance all orbits have approximately the same mean value for 
I\2- As we show below, this can be improved by increasing the 
number of fitted positions on orbits to improve the sampling of 
the various tori of orbits, and by increasing the order of the fit- 
ting polynomial to allow the modelling of tori with smaller scale 
structures. 



Semi-ergodic orbits are also present and are visible, confined 
between the tube and box orbits at E - -0.05 (Figure [5]). At en- 
ergies E = -0.05, they correspond to orbits that pass close to 
the centre with x^ from to ~0.07. Some of these orbits are 
included in the fit but do not alter the fitting quality: a possible 
cause is that they remain confined within a small volume and 
thus are reasonably adjusted as the immediately nearby regular 
orbits. The dispersion of /12 of the semi-ergodic orbits is signifi- 
cantly higher than for the nearby regular orbits, however. 

Figure [5] (top-right) shows the Poincare surface of section at 
E = -0.05 built from the fitted integral of motion /12. We note 
that only two families of orbits are identifiable because no high- 
order resonances are modelled. 



3.2. Higher orders 

To improve the approximate integral using the same data 
within the same time interval Ar = 4000, we progressively 
increase the order of the fitted polynomial up to order 22 (2057 
coefficients). By improving the fit, we mean that we reduce 
the variation of / within a given time interval Ar. In contrast, 
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we remark that formal forms of the second integral obtained 



also with polynomial series (Whittaker 1937 Birkhoff 1927 
Contopoulos|| 1 960 Gustavson 1966 1 lead to polynomial series 



that are asymptotic, i.e., up to an optimal order the accuracy 
of the series improves, and beyond this optimal order the 
accuracy decreases, the series being divergent. This is related 
to the Nekhoroshev| ( |1977[ ) theory (see also Contopoulos 2002 
§2.3.6), which establishes that the formal integrals are valid 
over exponentially long (or short) times. Here, our approximate 
integral is different in the sense that it is always adjusted for 
a finite time interval, and its validity for longer time intervals 
has to be checked (the accuracy within the finite time interval 
can always be improved up to a certain limit L, and we expect 
that this limit L will increase when increasing the time interval, 
so that it does not contradict with the Nekhoroshev theorem.) 
With the 2D potentials considered here, the long-term behaviour 
of the orbits can be constrained by absolute barriers in the 
phase space due to KAM tori (at the opposite of the case of 
higher dimensional potentials that allow Arnold diffusion). This 
may imply that our approximate integral variations would not 
increase for longer time. However, in the case of semi-ergodic 
orbits, orbits may be confined within a restricted domain of 
the phase space before reaching another domain if the domains 
are scarcely connected, as illustrated by Athanassoula et al. 
( |1983| figure 8). In such cases we would expect that the integral 
variations increase for longer times. 

Figure [6] shows the mean / and dispersion 07 for each orbit, 
the values are presented sequentially by group of energy (left 
to right E varies from -0.95 to -0.05). Increasing the polyno- 
mial order allows us to model more resonances. At order 10 the 
only one minimum of / corresponds to the resonant orbit 1:1. 
At higher polynomial order, / has more extrema correspond- 
ing to secondary resonances. The dispersion through these res- 
onances decreases progressively as the order of the polynomial 
increases: from order 12 to 20 the mode (histogram maximum) 
of 07 changes from 0.050 to 0.0007, while the mean of 07 varies 
from 0.18 to 0.011. The mean dispersion increases at order 22, 
but the highest values of 07 at resonances decrease significantly. 

We define a more quantitative criterion to consider the ability 
to distinguish two nearby orbits with the same energy: the rela- 
tive dispersion as the ratio of the dispersion 07 to the variation 
of / between orbits: 

CT re l = <T I /\dI/dx\ E . 

This relative dispersion has the dimension of x, here with 
values of a few 10~ 4 and can be compared to the full range 
of variations of x from to about 2. This relative accuracy is 
the poorest at resonances when (dI/dx) E -0, but the mode of 
cr re i varies from 5. 10~ 4 to 1.5 10~ 4 when the polynomial order 
changes from order 12 to 20. The mean cr re i (clipped within [0 
to 1. 10~ 2 ]) decreases first but then flattens at 1.0 10~ 3 from the 
orders 14 to 20. 

Increasing the polynomial order allows us to model more res- 
onances. At low order from 8 to 16 only the two main families 
are modelled, as is visible in the reconstructed Poincare section 
at E — -0.05 (Figure [5j. At orders 18 to 22 three more reso- 
nances are visible in the reconstructed Poincare surface of sec- 
tion corresponding to families 2:3, 3:4, and 3:5. Only the family 
2:3 (extended from x — 0, u — 1.06 to x = 1.9, u = 0) is ac- 
curately modelled, the two other ones are incorrectly positioned. 
Considering the small size of the structures of these resonances, 
we suspect that a higher order polynomial fit is necessary and 
also because increasing the polynomial order progressively im- 
proves the fit. 




1000 1500 
Sequence 




Fig. 6. Integral / for 2085 orbits covering 19 energies. From top 
to bottom, fits with polynomials of order 10,16 and 22. Orbits 
are plotted sequentially according to their initial position *,,„,, 
and according to the energy (E = -0.95 left to E = -0.05 right). 
The red line is the dispersion of / multiplied by 10. 



Owing to limited computing resources we did not ex- 
plore this systematically with much higher polynomial orders. 
Nevertheless, we explored fits with a restriction to orbits just 
close to the energy E = -0.05 using 2200 orbits. The resulting 
fit is better (Figure [5} with now two resonances correctly posi- 
tioned within the surface of section, but only the family 2:3 has 
the correct multiplicity of islands. Now the family 3:4 is also 
correctly positioned, but the number of reconstructed islands is 
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incorrect, but by increasing the polynomial order, the exceeding 
island progressively disappears. 

Finally, we proceed to the local fit of the three resonances. 
Figure [5] shows the reconstructed Poincare surface of section 
from the three analytical approximate integrals of motion that 
shows a satisfying fit. A second 1 :4 family (not used for the fit) is 
also correctly modelled (the corresponding periodic orbit never 
cross the x-axis perpendicularly). 

All these results point out that for this potential, an approx- 
imate integral is constant with a very high accuracy. There is 
also a strong indication that a better modelling of / could be 
achieved by using a polynomial of sufficiently high order. 



3.3. Other combination of coefficients 

We have tried about 30 different combinations of coefficients to 
numerically build quasi integrals /: all with the same number of 
coefficients, and the fixed term is either x 2 , v 2 or x 2 v 2 . The effi- 
ciency or accuracy of 7 is determined according to the value cr re i 
for every orbit. The various combinations are not equivalent in 
terms of efficiency, some give a better fit for orbits with high en- 
ergies, others for low energies. The four following combinations 
give the best results 

- coefficient of x 2 set to 1 , all others v 2n and x 2k are zero, 

- coefficient of x 2 set to 1, all others y 21 and x 2k are zero. 
With these combinations we have 

7(0, 0, 0, 0, v) = and I(X sup , E , 0, 0, 0) = X 2 supE 

- coefficient of v 2 set to 1 , all other v 2n and x 2k are zero, 

- coefficient of v 2 set to 1, all other v 2n and u 2m are zero. 
With these combinations 

7(0, 0, 0, v) = 2(E - O(0, 0)) and I(X sup , E , 0, 0, 0) = 0. 

Applying the permutations of coefficients x <-> y and also 
u <-> v gives satisfying adjustments. 



4. Logarithmic potentials 

To check that the numerical results described in the previous 
section have some generality, we considered other potentials. 
These potentials have a very limited amount of ergodicity, the 
chaotic orbits remaining confined within a small volume of the 
phase space. We also limited most of our fits using polynomials 
with order smaller than 18 to avoid prohibitively long computing 
times (100 h.cpu). 



4.1. Scale-free logarithmic potential 



The logarithmic potential introduced by Richstone (1980 1 to 
model axisymmetric disc galaxies produces a flat rotation curve: 

<& = log(x 2 +y 2 /q 2 ). 

The parameter q allows us to model a flattening of the poten- 
tial and of the corresponding density. The density is not positive 
everywhere when with q < V2/2, and that with q < V3/4 ~ 
0.87 the density isocontours present a depression close to the 
y-axis. 

A detailed analysis of orbits in this potential was pre- 
sented by Miralda-Escude & Schwarzschild (1989) with plot- 
ted Poincare surfaces of section. This potential has semi-ergodic 




Fig. 7. Logarithmic potential: (x, u) Poincare surface of section 
at E = for three families of orbits (circles) and three different 
analytical tori (continuous lines) obtained with polynomials of 
order 18. 



motions ( Papaphilippou & Laskar 1996) and is not integrable. 
Chaotic orbits have excursions close to the centre and circulate 
around the sets of tori of the major families of resonant orbits. 
The phase space is mainly occupied by regular orbits (at least 
with q from 0.2 to 1) and the dominant families of these regular 
orbits are tube orbits of type 1:1. Box orbits of type 2:1, 3:2, 4:3, 
etc (also named boxlets orbits : banana, fish, pretzel, etc) cover 
most of the remaining surface of the Poincare section. 

Using the procedure described in Section 2, we determined 
polynomial forms to evaluate a second integral of motion by fit- 
ting orbits of type 1:1, 2:1, and 3:2 (this is done only with E = 
since the potential is scale-free and the integral can be simply 
deduced for other energies). Using polynomials of order 18, we 
succeeded in modelling most of the orbits of each of these three 
families separately. Thus, a different polynomial was used for 
each of the three families (Figure [7]). We also obtained a gen- 
eral fit of the three families with a unique polynomial of order 
28 (4760 coefficients), but just two-thirds of the orbits plotted 
in Figure [7] are fitted (the orbits close the central periodic or- 
bits). Within the surface of section (x, u), these unfitted orbits 
are poorly fitted when x ^ 0. 1 but are correctly fitted at larger x. 

The scale-free logarithmic potential is the one for which the 
application of our method gives the most limited success. As 
for the other potentials, we used our program to fit the energy 
E(x,y, it . v) considering it as a polynomial and consequently also 
the potential <5(x, y). Not surprisingly, a polynomial form hardly 
fits a logarithmic function without core. Furthermore, examining 
the Boltzmann equation for this potential, we can assume that to 
build a second integral with a formal series, a series of ratio- 
nal functions will be more suited. Here, the Prendergast (1982) 
method to obtain rational solutions should be certainly more ef- 
ficient. 

4.2. Logarithmic potential with a core 



As suggested by Richstone (1980) and already used by Binney 
( |1981| l andpvlagnenat|( |1982| l, the logarithmic potential 



<E> = log(a 2 + x 2 + y 2 /q 2 ) - 2 log(o) 
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Fig. 8. Logarithmic potential with a core: (x, u) (top) and (y, v) (bottom) Poincare surfaces of section at £ = 3. Fits with an 18 
order polynomial from orbits with E — to 3.5. Computed orbits (right) and anaytical tori (left). 



has a core. We examine this potential with a — 1 and q = 0.8. 
The presence of a core drastically reduces the extent of the semi- 
ergodic orbits for energies E < 4 (and a = 1). It is at much 
higher energies and with x or y coordinates much larger than 
the core radius a, that orbits are similar to the orbits within the 
logarithmic potential seen in Section 4.1. At low energies box 
orbits dominate, type 1:1 tube orbits appear at E ~ 0.9. 

We obtain a good fit for orbits from E = to 3 with a I2' h or- 
der polynomial (and coefficients selected as in Section 3.1), the 
fit is good everywhere except close to the transition between box 
and tube orbits, the less accurately fitted region being at the ap- 
parition of tube orbits when E ~ 0.9. The fit is greatly improved, 
however, by increasing the number of orbits in this transition re- 
gion, for instance using 23428 orbits (1000 positions along each 
orbit) and order 16. Some resonant orbits are not modelled (see 
Figures [8] (y,u) Poincare section at E = 3.0, and (y, v)=(0.8,0) 
or (0,1.5) or (1.7,1.5)) but they cover a very small volume of the 
phase space. Increasing the order to the 18' /! order, the fit im- 
proves, but high-order resonances are still poorly fitted. Fitting 
orbits from E = to 3.5 (with 27619 orbits) with 5000 positions 
on each orbit and order 18, the fit improves from E — to 3. 

4.3. Axisymmetric scale-free logarithmic potential 

Still following Richstone (1980), we consider a 3D flattened 
logarithmic potential and motions with non-null angular mo- 
mentum. We set no core (a = 0) and a flattening parameter 
q = 0.8 and consider only orbits with the same angular momen- 



tum L, — I. Results for any other non-null L~ are easily deduced. 
Thus we consider orbits in the effective potential 

4>eff = f log(tf 2 + Z 2 /q 2 ) + ^2 - f (log L z + -). 

E e ff close to zero corresponds to nearly circular obits around 
the z-axis. The radial orbits with z — and a null vertical velocity 
dispersion <x, are a first family of periodic orbits. The shell orbits 
with a null radial velocity vr when z — are another family. 
For a given E e g, the amount of radial versus vertical extension 
of an orbit is constrained by a second integral (see for instance 
Ollongren 1962). The rotation prevents orbits from circulating 
close to the centre and thus reduces the amount of ergodic orbits 
by comparison to the case with L z = 0. Chaotic orbits occupy an 
extremely small volume in the phase space, at least for E e g < 8. 
At low energies (E 5= 1) orbits are mainly box orbits. 

Radial orbits have z = and at their two extrema R m i n 
and /? max have also vr = v z = 0. Thus, for all such orbit: 
^(^min> 0,0,0) = I(R m i„, 0,0,0), a relation that can be achieved 
only if I(R, 0,0,0) = f(®(R,z = 0)). The choice of coefficients 
in Section 3.1 does not allow us to fulfil this condition. Therefore 
here we set to one the coefficient of v? and remove all v 2 ", with 
n > 1, and all terms (removing all v?' 1 , with m > 1, and all 
R 2k is also satisfying.) 

Figure [9] shows the results of a fit with 24000 orbits (100 
positions on each orbit) with E — to 0.5, and an 18 f/l order 
polynomial. 
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Fig. 9. Axisymmetric logarithmic potential: orbits with L z — 1. 
(top to bottom): (x, u) Poincare surface of section with E = 0.01 
and .5. Computed orbits (dotted red circles) and analytical re- 
constructed tori (continuous dark lines) with an 18'' 1 order poly- 
nomial fitted along orbits with E = to 0.5. 



Fig. 10. Toda lattice: (Top) (y, v) Poincare section at (x-0, u > 0) 
and £=1/12 and analytical tori (continuous lines) with a \A th 
order polynomial. (Bottom) (y, v) Poincare section at (x-0, u>0) 
and E = 1 for a few orbits (dotted red circles) and analytical tori 
(continuous dark lines) with a 16 th order polynomial. 



5. Henon & Heiles and Toda lattice potentials 

We do not present the various tests of our procedure performed 
with integrable dynamical systems: axisymmetric, Stackel po- 
tentials, or systems with a second integral with a finite polyno- 
mial form. 

Instead, we summarise the results obtained for two classic 



potentials, the |Toda ( 1970| > lattice (2D case) and the |Henon 
|& Heiles] ( |1964| l potential (a historical description concerning 
these potentials on simulations in dynamics is given by |Weissert 
1997\ . The Toda lattice potential has a known analytical second 



integral ( |Henon|[T974l |Manakov|[T974l [T975] |Flaschka| [1974) 1 



The Henon & Heiles potential has regular orbits at low energies 
while nearly full ergodicity appears at high energies. The two 
potentials are closely linked since the first four terms of the 
series development of the 2D Toda lattice give the Henon & 
Heiles potential, one is integrable, the other one shows a nearly 
complete dynamical chaos at high energies. The functional 
form of the second integral of the Toda lattice is known and 
is odd for the momenta reflecting that there is no box orbit 
but only circulating motions (tube orbits). The extrema of the 



second integral at fixed energy allows the precise identification 
of the family of stable orbits. This potential and its second 
integral can be written as a polynomial series with an infinite 
radius of convergence for all coordinates. This was assumed 
to be important as a first step to test the procedure developed 
in this paper, which assumes that a polynomial series exists to 
represent the hypothetic second integral. 



5.1. Toda lattice 

Selecting 79 orbits with energy E = 1/12, 10000 positions on 
each orbit, and a fit at order 14, we obtain a median relative 
accuracy of 1.2 10~ 5 . We obtain 7.5 10~ 7 at order 16 (Fig. 10 1. 
Selecting a uniform distribution of 1700 orbits from energy 
E = to 1, still with 10000 positions, we obtain a median rel- 
ative accuracy of 1.5 10 _I at order 14 and 3.5 10~ 2 at order 18. 
We also succeed in obtaining a good fit from E = to 4 with 
2800 orbits. The polynomial expansion of the known Todda lat- 
tice second integral contains only terms of order 0, 1 or 3 in 
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Fig. 11. Henon & Heiles potential: (Top) (y, v) Poincare section 
at (x=0 and u > 0) and E - 1/12 for a few orbits (dotted red 
circles) and analytical tori (continuous dark lines) with a 12'' 1 or- 
der polynomial. (Bottom) Poincare section at E — 1/8 for a few 
regular orbits (dotted red circles) and analytical tori (continuous 
dark lines) with an 18'* order polynomial, the unplotted domain 
is covered by chaotic orbits. 



the momenta v. Our a priori construction and fit of a polyno- 
mial series does not impose such a constraint on the coefficients, 
accordingly our fitted series requires a much larger number of 
terms and coefficients to achieve the same accuracy as using the 
polynomial form of the known integral. 



5.2. Henon & Heiles potential 

Within the Henon & Heiles potential, we compute 10000 posi- 
tions along 41 orbits with the energy E = 1/12 (AT=10000). 
The recovered analytical tori are quantitatively correct and bet- 
ter than a few thousandths in relative accuracy at orders higher 
than the 1 th (83 coefficients). Figure 
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may be compared to the 
result of Gustavson (1966, figure 13), where the integral of mo- 
tion was determined by symbolically manipulating polynomials 
up to the order 7. At order 12, our median relative accuracy is 
improved: 2 10~ 4 , and at order 14 it is 10~ 4 . 



At energy E — 1/8 an ergodic orbit covers about half of the 
phase space (see figure 5 in |Henon & Heiles||1964| ). We obtain 
a roughly approximate integral of motion for the regular orbits 
close to the main periodic orbits. With a fit of 54 regular orbits 



close to the resonance 1:1, the recovered analytical tori (Fig. 1 1 
show that the three main families are recovered. This may be 
compared to the results of Gustavson ( 1966, figure 10) obtained 
at the same energy: he obtained analytical sections in good quali- 
tative agreement with the sections plotted from computed orbits. 



More recent works (see 


Kaluza & Robnik 1992; Robnik 1993; 


Contopoulos et al. ['2003 


i give higher orders of formal integrals 



and show a better agreement at energy E=l/8 where the chaotic 
orbits occupy a large volume of the phase space. 

6. Conclusion 



Obtaining integrals of motion of a dynamical system simpli- 
fies the analysis and the understanding of the system, because 
it allows us for instance to build more easily tractable distribu- 
tion functions. It also gives a synthetic description of the orbits, 
the building blocks of galaxies. This is well-known for axisym- 
metric systems, and it has been also extensively developed for 
Stackel potentials. Other systems with known analytic integrals 
are apparently less interesting for the astronomical community. 

Therefore many efforts and methods have been developed to 
obtain approximate integrals of motions for more general poten- 
tials. Moser's theorem states that most invariant curves will be 
preserved under a weak perturbation of an integrable dynami- 
cal system. A method initiated by Birkhoffj ( |1927[ ) consists of 
writing the Hamiltonian in the so-called normal form, which al- 
lows the formal construction of polynomial series that are so- 
lutions of the Boltzmann equation.These series are generally di- 
vergent (Siege! 1941 Contopoulos et al. 2003 [l. Even for an inte - 
grable system, they may be not convergent ( Wood & Ali| 1987 >. 



Gustavson ( 1966) developed an algorithm to obtain these formal 



integrals and applied it to the Henon & Heyles potential. 

A more direct method by Whittaker ( |1937| l or |Contopoulos| 
( 1960) also consisted of looking for formal polynomial integrals 
of motion, but did this directly by substituting and comparing 
coefficients. Giorgilli & Gal gani| ( |1978) > solved the consistency 
problem of these direct methods. 



Inspired by these methods, Giorgilli & Galgani ( 1978 1 pro 



posed a new method for constructing formal integrals near an 
equilibrium point, and a numerical program (Giorgilli 1979]> was 
used for applications (see for instance Kaluza & Robnik 1992 



[Contopoulos et al. 2003). The formal integrals are different for 
resonant and each non-resonant orbit ( |Contopoulos et aLp OOO). 

A different technique, called torus construction (for a review, 
see |Valluri & M erritt 1999 ), consists of numerically mapping the 
action-angle coordinates of a known potential into the action- 
angle coordinates for the system under investigation (McGill & 
|Binney|1990| ). This is obtained by a least-squares procedure. It 
implies the availability of a close and integrable Hamiltonian 
whose tori can be simply mapped to the target tori; it also im- 
plies mappings from different toy tori for different orbit fami- 
lies (Kaasalain en & Binney|19 94). Moreover, there are general 
methods of semi-numerical perturbation that take into account 
the full 'distortion' of the invariant tori ( |Henrard[l 990). 

The method developed in this paper partly builds on pre- 
viously published methods. We solved the Boltzmann equation 
with a polynomial, but numerically and with a fit to many pecu- 
liar solutions (i.e. orbits). This method is a priori limited by the 
accuracy with which the invariant curves are defined on the sur- 
faces of section, i.e. to which point they are effectively invariant 
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and affected by diffusion. This limiting accuracy depends, for a 
given potential and a given orbit, on the time interval of integra- 
tion of the orbit. But for a given time interval, it should be suf- 
ficient to increase the order of the fitting polynomial to achieve 
higher accuracy. We succeeded in different examples to fit vari- 
ous resonant and non-resonant families of orbits with the same 
integral. The quality of the fits depends on selected coefficients 
of a polynomial series. Here, the method was applied to various 
2D potentials representative of motions within elliptical galaxies 
or motions within a non-axisymmetric disc. For all these poten- 
tials, a unique analytical integral of motion was obtained for an 
extended range of energies and regular orbits. The coefficients 
of the polynomial of the integrals were obtained from a linear 



least-squares minimisation. |Warnock ( 1991 1 and Bazzani et al. 



( 1991 ) also developed methods of orbit fitting but with limita- 
tions; the former with restrictions to high-order resonant orbits 
or non-resonant orbits, the latter to 2D sections. In addition, they 
did not impose the energy independence of their integral, which 
may impact the numerical accuracy of the fitted second integral. 

As a future development, the method presented in this pa- 
per will be applied to the construction of a simple 3D galactic 
model. We will use a scale-free logarithmic disc-halo potential 
(Bienayme 2009 1, for which a third integral h will be deter- 
mined: the third integral is first determined for a given value of 
the energy and a given non-zero angular momentum, and is then 
simply deduced for other values of E and L z . Afterwards, the 
action can be evaluated from I3. Finally, a distribution function 
for an exponential stellar disc can be built following the general- 
isation of the Shu distribution function by |Kuijken & Dub inski 



( 1995|l, [Bienayme & Sechaud]( [T997l l, or |Bienayme|dl999| ). This 



model will allow us to explore various properties: for instance, 
to check the domain of validity of the asymmetric drift relation 
especially out of the mid-plane, to examine the velocity ellipsoid 
tilt dependence on the distribution function, or to build a realis- 
tic model to measure the K z force from data far above and below 
the galactic plane by minimising the number of free parameters 
and assumptions in the modelling. 
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